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Abstract 

The sparsity and algebraic structure of convex resource allocation problems have given rise 
to highly specialized methods. We show that the sparsity structure alone yields a closed-form 
Newton step for the generic primal-dual interior point method. Computational tests indicate 
that this makes the interior point method competitive with the leading special-purpose method 
on several important problem classes for which the latter can exploit additional algebraic struc- 
ture. Moreover, the interior point method consistently outperforms the specialized method when 
no additional algebraic structure is used or available. 

Keywords: convex programming; interior point methods; continuous knapsack 

1 Introduction 

The continuous nonlinear resource allocation problem (or continuous nonlinear knapsack problem) 
is a finite-dimensional convex optimization problem involving a separable objective, simple bounds 
on each variable, and a single constraint linking the variables by means of a separable function. 
The survey paper of Patriksson [25] shows that the resource allocation problem has a long history 
and diverse applications; see also the book by Ibaraki and Katoh [14]. The contexts in which the 
problem appears often demand that it be solved very quickly, even when the problem dimension 
is quite large. Consequently, researchers long ago moved beyond general-purpose nonlinear or 
convex programming procedures and focused on exploiting the special structure of the optimality 
conditions for the problem. 

As noted by Patriksson [25], two frameworks have emerged as the most competitive for solv- 
ing resource allocation problems: the pegging or variable- fixing methods and the breakpoint-search 
methods. Each involves solving a finite sequence of optimization subproblems. In a pegging method, 
the subproblems correspond to ignoring the bounds on some variables while holding other variables 
fixed ("pegged") at one of their bounds. Each iteration identifies at least one variable that can 
be pegged, and it stays pegged thereafter. The procedure ends when all unpegged variables (if 
any remain) in some subproblem satisfy their bounds at optimality when those bounds are not 
explicitly enforced. By contrast, a breakpoint search acts within the one-dimensional dual space 
corresponding to the Lagrange multiplier for the single linking constraint. The concave dual objec- 
tive can be viewed as a piecewise-defined function whose breakpoints are easily listed, so a binary 
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search identifies either a maximizer among those breakpoints or a pair of breakpoints that most 
closely bracket a maximizer. At each tested breakpoint, a separable optimization is performed 
in the primal variables; the final interpolating subproblem is essentially equivalent to a pegging 
subproblem. 

The pegging and breakpoint-search frameworks have each been rediscovered many times, so 
original attribution is difficult. Patriksson [25] cites Theil and van de Panne [33] as a possible 
genesis for the pegging methods, with the first formal description apparently given by Sanathanan 
[29j . A dual approach seems to have been described first by Churchman et al. [8], with the earliest 
explicit search of breakpoints carried out by Charnes and Cooper jTJ. In summarizing the state of 
the art, Patriksson [25] notes that computational studies in the literature have generally indicated 
that pegging is superior to breakpoint search when certain subproblems (see the next section) 
common to both methods are easily solved, whereas breakpoint search is faster otherwise. He also 
observes that numerical comparisons of either method with general-purpose solvers are essentially 
nonexistent in the literature. 

In the current paper we present computational results showing that a primal-dual interior point 
method outperforms breakpoint search on problems where the latter is traditionally considered to be 
the best possible choice, namely, when its subproblems do not admit closed-form solutions and must 
be solved numerically. Each iteration of a primal-dual interior point method first calculates a search 
direction for Newton's method applied to a perturbation of the Karush-Kuhn- Tucker conditions, 
and then chooses a step-size so that the primal and dual variables satisfy their respective bounds 
strictly [23]. The key to the method's success in the present setting is the special form of the linear 
system that is solved when calculating the search direction. Aside from exploiting the linear-system 
structure, our implementation of the interior point method has been kept fairly simple and generic. 
By contrast, we have taken pains to make the breakpoint search as efficient and reliable as possible, 
so as to justify our conclusions about the superiority of interior point methods in this setting. 

In the next section, we state the problem formally and review its optimality conditions. In §0 
we describe the breakpoint search and interior point methods, along with details of their implemen- 
tation. Section [4] lays out the problem instances used for the computational tests, and the results 
are discussed in £}5j Some concluding remarks are provided in £j6) 

2 Notation, assumptions, and optimality conditions 

We consider the resource allocation problem in the form 

n 

(1) minimize f(x) := ^fijxj) 

i=l 
n 

(2) subject to g(x) := ^ gi(xi) 

i=i 

(3) I < x <u. 

Here x, I, and u are n-vectors of real numbers, whereas b is a real scalar. Inequalities of vectors 
are interpreted coordinate- wise. For the purpose of the present study, we make the following 
assumptions: 

Al. The functions fi and gi are convex and twice differentiable on an open set containing the 
interval [k,Ui\. 

A2. The relaxed problem, in which ([2]) is replaced by g(x) < b, has no optimum with g(x) < b. 



over all x 
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A3. The function /j is decreasing on and gi is increasing on with g(l) < b < g(u). 

The randomly generated test instances discussed in later sections all satisfy these assumptions. 

In practice, we are actually interested in the relaxed problem of Assumption A2. However, 
we can cheaply determine whether that assumption holds provided we know the intervals of mono- 
tonicity for each fa. Indeed, most presentations of methods for resource allocation problems include 
Assumption A2 because it can be enforced through some combination of initialization, preprocess- 
ing, and data generation. The required computational time is typically quite small compared to the 
time needed to solve the problem; in many instances, the task is trivial. Likewise, many treatments 
also include some form of Assumption A3: it simplifies discussion and implementation, and it can 
be guaranteed inexpensively by first identifying the intervals of monotonicity for /j and gi and then, 
as needed, shrinking the intervals and reorienting the variables Xj. 

Assumptions A1-A3 imply that problem ([l])-([3]) and the relaxed problem are equivalent and 
admit an optimal solution, and they also guarantee that the Slater constraint qualification holds 
for the relaxed problem. By Lagrangian duality, necessary and sufficient optimality conditions for 
([I])-© can therefore be expressed as follows: g(x) = b and, for some real number p, x is a solution 
to the separable optimization subproblem 

(4) minimize f{x) + pg{x) over all x subject to ([3]). 
The dual objective is 

n 

(5) p^-bp + Y] mm [fi(xi) + pgi(xi)}, 

• Xie [k,ui] 

which attains its maximum; moreover, any maximizer p is necessarily nonnegative. The subproblem 
(j3J has coordinate-wise optimality conditions given by 

f'i (xi) + pg'i(xi) = 0, if k < Xi < Ui, 
f-(xi) + pg'i(xi) > 0, if Xi = l h 
fi(xi) + pg'i(xi) < 0, if xi = Ui. 

The left-hand sides give the Karush-Kuhn- Tucker multipliers for the bounds k < X{ and Xi < Ui, 
respectively, as 

Xi := max{0, ~[fU x i) + P9i( x i)]} and Mi : = max{0, f-(xi) + pg'^Xi)}. 

Letting s := u — x denote the vector of slack variables for the upper bounds on x, we can express 
the Karush-Kuhn- Tucker (KKT) conditions for ([T])-([3]) as 

(6) V/(:c) + pVg(x) - A + n = 0, 

(7) x + s = u, 

(8) x > I, A > 0, s > 0, p > 0, 

(9) diag(x — l)X = 0, diag(s)/i = 0, 

(10) g(x) = b. 

Here the notation diag(z) refers to the diagonal matrix whose diagonal entries are given by the 
entries in the vector z. An equivalent version of the Lagrangian dual incorporates the multiplier 
vectors A and p, but that form is not used in the sequel. 

The three solution frameworks discussed in the introduction exploit the above developments in 
different ways: 
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• Pegging methods solve subproblems of the form (H])-©, but for which some variables are 
held fixed while the bounds ([3]) for all remaining variables are omitted. 

• Breakpoint search methods maximize the dual objective ([5]) by solving a sequence of sub- 
problems of the form @ at various values of p. 

• Primal-dual interior point methods apply Newton's method to perturbations of the KKT 
system P- (fT()]) . 

The pegging and breakpoint search methods both benefit considerably when minimization of Xi \— > 
fi{xi) + pgi(xi) can be handled efficiently. As noted by Patriksson [25], previous computational 
studies have indicated that pegging is generally the fastest method for resource allocation problems 
in situations where the pegging subproblem can be solved in closed form, and that breakpoint search 
is the best method otherwise. The present article shows that, in fact, interior point methods are 
competitive with breakpoint search in the latter situation. Because we focus on problems for which 
breakpoint search dominates pegging, we do not include pegging methods in this study. Indeed, 
the pegging approach is not even well-defined for some of the problems we consider, because the 
pegging subproblems do not admit optimal solutions. 

3 The methods and their implementation 

In this section, we describe the two main approaches tested in our computational study. In addition 
to the details indicated below, our MATLAB implementations of both methods employ vectorized 
operations with logical indexing wherever they provide an efficiency gain over explicit loops or 
conditionals. 

3.1 Breakpoint search 

Breakpoint search (see Figure [1]) is based on the observation that the dual objective ([5]) is concave 
and defined piecewise with a finite number of easily calculated breakpoints. The derivative, or 
subdifferential, of this objective is nonincreasing. A binary search of the breakpoints therefore 
identifies either one that is itself a root or a pair of breakpoints that most closely bracket a root. 

There are at most 2n breakpoints, occurring at p- values where some Xi h-> fi(xi)+ pgi(xi) attains 
its minimum over [Zj,itj] at an endpoint l\ or U{. Equivalently, a breakpoint makes the derivative 
x i ^ fl( x i) + P9i( x i) nonnegative at ij or nonpositive at u\. Consequently, all breakpoints have the 
form 



The monotonicity of fi and §i allow us to define pf = oo when g[(li) = and to guarantee that 
g[{ui) > in the definition of pj . The convexity and monotonicity of fi and gi also guarantee that 



The binary search sequentially refines a bracketing p~ < p* < p + until the true root p* lies 
between two consecutive breakpoints. The bracket is adjusted inward by extracting a breakpoint p 
within it and testing the sign of the derivative of the dual objective ([5]). To evaluate that derivative 
at p, we first fix 




fl(k)/g'i(k) or Pi := -f!(ui)/gfi(ui). 



< PT < Pj- 



(11) 
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Input: data /, g, b, I, u defining an instance of problem (HJ) — ([3j) satisfying assump- 
tions A1-A3 of g2j 

Output: approximate solution x for problem (d])-©. 

1. [Initialization] 

Set I := {1, . . . , n}, p~ := 0, p + := oo and b := b. 

For all i G I, set p+ : = and P'i : = "4rT- 

Set i? := {Pj 7 }ie/ U {/^~}j g /, stored as an unsorted list of possibly replicated 
values. 

2. [Breakpoint selection] 

Select p G i? as a median value from the list representing i?. 

3. [Subproblem solution] 

Set I~ := {i G / : p < p^} and /3~ : = ^ 9i(ui). 

iel- 

Set /+ := {i G / : p > pf} and /3+ := ^ 

iej+ 

For each i £ I \ (I~ U find Xj to minimize fi(xi) + pgi(xi) subject to 

JiA < ~~^ • 

Set p := p~ + 0+ + ^(sf). 

iei\(l-ui+) 

4. [Bracket update] 
Set Ab := 0. 

If /3 > b, then: 

Set p" := p, Ab := /3~, I := I \ I~ and R := {r G R : r > p}. 
For each i G I~, set Xi := Uj. 

If p < b, then: 

Set p+ := p, A6 := Ab + 0+, I := / \ /+ and R := {r G : r < p}. 
For each i G I + , set Xj := l; L . 

Set S := b - Ab. 

5. [Stopping criterion] 

If R 7^ 0, then go to stepEJ Otherwise, continue to step El 

6. [Interpolation within final bracket] 

Find p G [p~ , p + ] and Xi for each i £ I so that fl(xi) + pg\(x,i) = and 



Figure 1: Breakpoint search 



The remaining minimizers are critical points: f[{xi) + pg\{xi) = and U < x% < ui. Depending on 
the problem data (see JQ), these critical points might be found (a) in closed form, (b) by using a 
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problem-specific implementation of Newton's method, or (c) by means of a general-purpose New- 
ton's method with Armijo linesearch for sufficient decrease and damping (as needed) to maintain 
li < Xi < U{. The derivative value at p is then given by — b + ^ i gi(xi), the sign of which determines 
whether p becomes the new p~ or p + . This in turn determines, through (jlip . that some values of 
Xi shall remain fixed and can therefore be removed from further consideration. 

The final bracket, if nontrivial, consists of two closest breakpoints with the optimal value of p 
lying somewhere between them. To interpolate between them, our implementation finds p and the 
unfixed Xj-coordinates (denoted by i G /) simultaneously by applying a multi-dimensional Newton's 
method with Armijo linesearch to the corresponding Lagrange multiplier conditions Y2i£i 9i( x i) = b 
and fl(xi) + pg'i(xi) = for i £ I. 

Throughout the procedure, the subproblem optimizations are initialized using the corresponding 
solutions from prior iterations. Also, we extract the required median values without sorting the list 
of breakpoints in advance, which can yield significant computational savings if each subproblem 
solution requires only a few operations per index i [6l [231 EH CE] • 



3.2 Interior point method 

The primal-dual interior point method solves the KKT optimality conditions (f6])- ([10p for the vari- 
ables (x, A, s, p, p). Its operation preserves strict inequality for the simple bounds ([8]), only allowing 
them to become active in the limit. The method is based on the perturbed KKT system 

(12) Vf(x) + pVg(x) - A + p = 0, 

(13) x + s = u, 

(14) diag(x — l)X = re, diag(s)/i = re, 

(15) g(x) = b, 

where e denotes the n-vector of all ones and the inequalities x > I, A>0, s > 0, p > are enforced 
separately. The system (fT2l - (fl~5j) is algebraically equivalent to the Lagrange multiplier equations 
for a related optimization problem involving barrier functions for the bounds ([3]): 

n 

minimize f{x) — r ^^[ln^ — li) + ln(sj)] over all x > I, s > 

i=l 

subject to g(x) = b, x + s = u. 

As the barrier parameter r > is driven to zero, we expect the (unique) solution (x, A, s, p, p) of 
(|12 p -(|15 p to tend toward the solution set of the original KKT system ([6])- (|10p , 

In each iteration of the interior point method, we calculate a Newton search direction for the 
perturbed system (fT2]) - (fT5|) and then take a step along that direction, damped so as to preserve 
x > I, A > 0, s > 0, p > 0. Next, the value of r is adjusted and the iteration repeats. The method 
is formally stated in Figure [2j using notation defined below. It is essentially a general-purpose 
interior point method for a nonlinear programming problem formulated with simple bounds and 
equality constraints. The notation emphasizes the single equality constraint of the present setting 
and the initialization exploits that aspect as well. Procedure parameters are given specific values 
(10~ 10 , 0.25, and 0.8 in steps 2, 3, and 5, respectively) that gave reliable performance in preliminary 
testing. 

The interior point method presented here is rudimentary and could likely be refined to provide 
theoretical guarantees of global convergence, superlinear convergence, or polynomiality in the num- 
ber of calls to /, g and their derivatives; see [SI EH EOl EH [HI dQl [32l EH ESI [35] for treatments 
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Input: data /, g, b, I, u defining an instance of problem (HJ) — ([3j) satisfying assump- 
tions A1-A3 of g2j 

Output: approximate solution x for problem (d])-©. 
1. [Initialization] 

Choose t G (0, 1) with g(tl + (1 — t)u) = b and set x := tl + (1 — t)u and 
s := u — x. 

Set p := 1, Aj := max{0, -[f< (x^+pg'^Xi)]} and Pi := max{0, fl(x i )+pg , i (x i )}. 

1. [Termination criterion] 

Set Td := Vf(x) + pVg(x) — A + p,, ri := diag(x — l)X, r u := diag(s)/i, and 

r g := g(x) - b. 

If none of the relative errors 

IMli 

1 + ||V/(a:)||i + \p\ + WVgWh + ||A||i + IMli ' 

IK||i \rg\ 

1 + \\x - + ||A||i' 1 + ||s||i + ||/i||i' 1 + \\g(x)\\x + \b\ 

exceeds 10 -10 , then we are done: x is considered optimal to within numerical 
tolerance. 

3. [Barrier parameter selection] 

Set a := — — and r := 0.25cr. 

In 

4. [Search direction] 

Set ri := n — i~e and r u := r u — re. 

Calculate (Ax, AA, As, Ap, Ap) to solve the linear system (|16|) . 

5. [Step size] 

Set t > to be the largest value for which 

(x,X,s,p) - i(Ax, AA, As, Ap) > (1, 0,0,0). 

Set t := min{l,0.8t}. 

6. [Update] 

Set (x, A, s, p, p) := (x, A, s, p, p) — t(Ax, AA, As, Ap, Ap). 
Go to step [2j 



Figure 2: Interior point method 



of such matters in similar contexts. However, the version presented here correctly solves all the 
instances described in $3] and handily outperforms breakpoint search on challenging problems of 
moderate to very large dimension. It therefore meets the needs of the present study. The key to 
making it competitive is the fact that the linear system in step U] can be solved in Cn arithmetic 
operations for a small fixed value of C, as we show next. 

To simplify the notation, we introduce a vector h with entries hi := f['(xi) + pg'/(xi) and let 
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£ denote x — I. We also use uppercase letters to denote these diagonal matrices: S := diag(£), 
A := diag(A), S := diag(s), M := diag(p), 77 := diag(/t). The linear system satisfied by the search 
direction is then 



(16) 
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To solve (fT6j) . first calculate the vectors 

w := h + H -1 A + S~V> 

z ■= W~ l Vg 

and let r\ := — lfVg T z. Notice that multiplying (I16p from the left by 



yields 
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From this we can read off the solution to (|16p as follows: 

Ap = 7?(r 5 - z T y), 
Ax = W x y - {Ap)z, 
As = -Ax, 
AA = H _1 n - E _1 AAx, 
A/i = 5 _1 r u - 5 _1 MAs. 

The values H _1 r/, H _1 A, S r u , 5 _1 M can be saved during the prior calculation of w and y 
and reused here. The solution of f)16|) therefore requires 9n — 1 additions/subtractions, 5n + 1 
multiplications, and 6n + 1 divisions. 

4 Test problems 

We classify the test problems considered in this work according to the availability of problem-specific 
subproblem solvers for breakpoint search. Namely, we place the problems into three categories as 
follows: 
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(a) Closed-form solutions are available for the separable breakpoint subproblem, but not for the 
breakpoint-interpolating (or pegging) subproblems. 

(b) No closed-form subproblem solutions are available, but each separable breakpoint subproblem 
has a readily identified initial guess that guarantees convergence of the standard Newton 
iteration. Typically, this guarantee is based on the knowledge of the convexity and inflection 
points of the nondecreasing function X{ i-> fl(xi) + pg'^Xi). 

(c) Neither of the above seem to apply. 

As described in the preceding section, we implemented the breakpoint search with a generic iterative 
solver for the subproblems. This was applied to all problem instances, regardless of category. 
However, we also tested implementations that exploit the problem-specific opportunities above 
when available. Specifically, for the first two categories above, breakpoint search was tested both 
with and without the problem-specific subproblem solvers; we therefore considered three approaches 
(including interior point methods) for those instances. For the third category above, only the two 
more general approaches were considered. 

Problems within each category are described in the subsections below. Instances were generated 
so that assumptions A1-A3 of $2] were satisfied, after possible reorientation of intervals. The 
notation z ~ U(a, b) indicates that the value z was selected according to a continuous uniform 
distribution on the open interval (a, 6), whereas z ~ N(/jl, a) indicates that z was selected according 
to a normal distribution with mean \x and standard deviation a. 

Some of the problem classes considered here have already been studied within the resource allo- 
cation literature, while others represent mathematical forms that commonly appear in operations 
research and its applications. We also draw attention to several important classes that have been 
left out: the quadratic knapsack P2 [12 El [Ml EEl 03] , stratified sampling [22l EOJ El [3] , and man- 
ufacturing capacity planning [HI 1361 El El E] problems. These admit closed- form solutions to the 
pegging subproblems described in £J21 and so the pegging method handily outperforms any of the 
methods investigated here. However, we mention in passing that the relative performance of our 
non-pegging methods on these three classes is very similar to what we report in ^5] for the nonlinear 
lot-sizing problem of §4.1.11 

4.1 Instances with closed-form solutions for breakpoint subproblems 

The problems in this subsection admit closed- form solutions to the breakpoint subproblem (J1J), but 
not to the interpolation or pegging subproblem. 

4.1.1 Nonlinear lot-sizing 

Problem in this class have fi(xi) = diXi + Cj/xj and gi(xi) = di/xi with domains X{ > 0, discussed 
in the book of Churchman, Ackoff and Arnoff |8j. Instances for the present study were generated 
as follows: 

• Gtj, Cj ~ U(l, 5) and a\ ~ U (1, 11); 

• k = \Jci/ai and m ~ U (k,li + 4); 

• b~U(g(l),g(u)). 
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4.1.2 Sum-of-powers objective with unit simplex constraint 

Problems in this class have fi(xi) = ai\Xi — yi\ Pi and gi(xi) = x«. Note that fi is everywhere twice 
differentiable when pi > 2. Instances were generated as follows: 

• ai~J7(l,10) and ~ [7(2,4); 

• k ~ C/(0, 5), nj ~ [7(Zj,Zj + 5) and ?/j ~ U(m,Ui + 5); 

• b~U(g(l),g(u)). 

4.1.3 Weighted p-norm objective with unit p-norm constraint 

Problems in this class have fi(xi) = ai\xi — yif and gi{x,i) = \xi\ p . Note /j and gi are everywhere 
twice differentiable when p > 2. Instances with p 6 {2,2.5,3,4} were generated as follows: 

• a* ~ 17(1, 10); 

• Zj ~ [7(0,5), iij ~ Z7(Zj,Zj + 5) and m ~ U(ui,Ui + 5); 

• b~U(g(!),g(u)). 

4.1.4 Decentralized inventory management 

This is the operational model of Bitran and Mondschein pQ. The constraint uses gi(xi) = CiXi and 
the objective function is given by 



where hi is a probability density on the real line. In our numerical tests, we take hi to be a normal 
density with mean m and standard deviation crj. Instances were generated as follows: 

• a = 0.2, (3 = 0.1, a = 0.35, and ck ~ 17(20, 120); 

• m ~ [7(5000, 15000) and ^ ~ [7(0.05^,0.10^); 

• Ui= min(argmin/j, with & ~ [7 (^ - 2<7i, 1.1 • argmin/j); 

• Zj ~ U(jm - 2ai,Q.l(jn - 2(n) + 0.9ui); 

• b~U(g(l),g(u)). 

Because f[ has bounded range, the Lagrange multiplier equation f[(xi) + pci = for the pegging 
subproblem often admits no solution. 

4.1.5 Reliability investment 

Problems in this class have fi(xi) = — ln(l — r^ 1 ) with domain x; t > and gi(xi) = QXi, as studied 
by Kettelle [15] and Everett |11| . Instances for the present study were generated as follows: 

• n ~ [7(0, 1) and a ~ [7(1, 10); 

• Zj ~ [7(0, 1) and ~ f7(Z;, 10); 
. b~U(g(l),g(u)). 
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4.2 Instances with good starting guesses for separable subproblems 

These problems have no closed- form subproblem solutions, but their structure provides a good 
starting guess for applying Newton's method to some reformulation of the breakpoint subproblem 



4.2.1 Resource renewal for processes with exponentially decaying throughput 

Problems in this class have fi(xi) = (iiXi{e~ 1 l Xi — 1) and gi{xi) = CiXi for X{ > 0, as studied by 
Melman and Rabinowitz [18J. For convenience, we extend fi to a C°° convex function on the real 
line by defining fi(xi) = —Xi for Xi < 0. Instances were generated as follows: 

• ai, Ci ~ C/(0.001, 1000) and take 7 = min^aj /cj}; 

{0, if ai/ci > 7, 

argmm/i(xi) + jgi{Xi), if cii/Ci < 7; 
Xi 

• Zj = and tij = 6/cj. 

The subproblem Q is solved efficiently by taking Xj = when 1 — pci /ai < and otherwise taking 
Xj = 1/yi, where yi > is obtained by applying Newton's method to the equation (1 + yi)e~ y — 
1 + pci/a>i = with initial value yi = 1. 

4.2.2 Weighted p-norm objective with unit r-norm constraint 

Problems in this class have fi(xi) = a.- L \xi — yi\ p and gi(xi) = \xi\ r . Note that fi and gi are 
everywhere twice differentiable when p, r > 2. Instances with p,r G {2,2.5,3,4} and p ^ r were 
generated as follows: 

• Oj ~ [/(1,10); 

• Zj ~ [7(0,5), Ui ~ [7(Zj,Zj + 5) and yi ~ U(ui,Ui + 5); 

• b~U(g(l),g(u)). 

The subproblem (j4]) is solved efficiently by applying Newton's method to //(xj) + pg[{xi) = 0, 
initialized at Xi = yi. 

4.2.3 Sum-of-powers objective with sum-of-powers constraint 

Problems in this class have /j(xj) = aj|xj — yi\ Pi and gi(xi) = \xi\ Vi . Instances were generated as 
follows: 

• Oi ~ [7(1,10) andpi,^ ~ [7(2,4); 

• Zj ~ [7(0,5), Uj ~ [7(Zj,Zj + 5) and yi ~ U(ui,Ui + 5); 

• b~U(g(l),g(u)). 

The subproblem (j4]) is solved efficiently by applying Newton's method to fl(x{) + pg\(xi) = 0, 
initialized at Xi = y^. 
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4.3 No special subproblem solution available 

The following problems don't appear to have a simple, specialized approach for solving the break- 
point subproblems. 

4.3.1 Convex quartic objective with unit simplex constraint 

This class of problems has fi{xi) = aixj + bixf + Cixf + dixi and gi(xi) = X{. Instances of these 
problems were generated as follows: 



• tH = (tf + vf)/V8, h = (C s -Ci + mXi)/V3 and a = {(? + xl)/V8, with fam, (i, Xi ~ N(0, 1); 
. d i = -f i {r i ), with r;~ £7(0,10); 

• Ui = min(Tj,Ai), with A, ~ t/(0,Tj); 

• k ~ f/(0,iti); 

• 6~?7( 5 (/), 5 ( U )). 



The choice of coefficients for /j guarantees that /j is strictly convex on the real line, which is true 
if and only if 8ajCj > 3bj, en > 0, and Cj > qO. Equivalently, the matrix 



must be positive definite. This can be ensured by selecting aj,6j,Q to be the rescaled entries of a 
matrix formed as A T A, where the entries of A are given by £j, rji, d, Xi- The choice of di guarantees 
that the critical point of /j is positive, after which the bounds are chosen so that each /j has the 
same monotonicity. 

It can be shown that the subproblem @ can be solved for these instances by applying Newton's 
method to fl(x{) + pg'^Xi) = 0, initialized at the inflection point cc, = — 6j/(4aj) of f[. However, 
such an approach turns out to be slower than using the general-purpose solver, which is why we 
have listed this class of problems here. 

4.3.2 Log-exponential objective with linear constraint 

Problems in this class have gi(xi) = CiXi and 



V8a,i V3bi 
\/3bi V8ci 
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Instances were generated as follows: 



• dij ~ N(0, 1) and a ~ 17(0, 10); 





if Zij > for all j, 
otherwise; 
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I min(argmin fi, 1.2Q ■ argmin /A if i 6 I, 

• k = Ui — 0.05|tij| — 5\rji\, with r/j ~ iV(0, 1); 

• b~U(g(l),g(u)). 

The rules used to generate aij, Ui and k ensure that fi is decreasing on [k,Ui], while also allowing for 
Ui to be a critical point of fi occasionally. In these instances, /■ might have multiple inflection points, 
so initialization of Newton's method for the breakpoint subproblem @ is not straightforward. 



5 Computational results 

The procedures of 33 were coded in MATLAB 7.12 and their performance compared on randomly 
generated instances of the problems described in 21 As noted, most of those problem classes admit 
some additional algebraic structure that potentially improves the performance of breakpoint search. 
Therefore, the results below refer to three optimization frameworks: interior point method, general 
breakpoint search, and specialized breakpoint search. All tests were performed on a dedicated 
2 x quad-core Intel 64-bit (2.26 GHz) platform with 24GB RAM running CentOS Linux. 

We attempted to solve instances of dimension 10 fc for k £ {1,2,3,4,5,6} with each of the 
applicable methods for each problem class. One hundred random instances were generated at 
each dimension for eight of the ten problem classes. The other two classes are those described 
in §4.1.31 and §4.2.21 for which we generated 100 instances at each dimension for each value of 
p,r £ {2,2.5,3,4}. Performance differences among these combinations of p and r were detectable, 
but too small to warrant separating out the results for discussion. We therefore report them in 
aggregate over all p and r, and remark that the higher values of p or r tended to require slightly 
longer running times than did the smaller values. 

Based on preliminary testing, an a priori time limit of 10 fc_2 seconds was imposed on each 
attempt at solution. The interior point method was run first and never exceeded this time limit. 
Consequently, for higher dimensional instances on some problem types, the breakpoint searches 
were limited to at most a factor of ten over the worst runtime for the interior point method on 
problems of the same type and size. 

The results of the tests are presented for each problem class in Figure [3) Because the breakpoint 
searches often exceeded the given time limits, the graphs consist mainly of the median runtimes. 
Mean runtimes are drawn separately whenever they can be visually distinguished from the medians 
on the scale shown; the mean curves are always the upper branches when two curves of the same 
line type are shown. Moreover, means that include runtimes at their limits are specially marked. 

Three notable features stand out in the graphs of Figure El First, the interior point method 
dominates the general breakpoint search whenever the dimension is at least 100. Second, the 
specialized breakpoint search dominates the general version by roughly an order of magnitude in 
many cases, but certainly not all (especially at higher dimensions). Third, the advantage of the 
specialized search over the interior point method tends to decrease with dimension, so that the 
interior point method emerges as dominant for high-dimensional instances of all but one of the ten 
problem classes. Table Q] displays the frequency with which such dominance occurs. 

Additional information regarding the distribution of solution times for the interior point method 
is given in Figured! Comparing that figure with Figure [31 we see that the variability among running 
times for the interior point method is less than the typical difference between it and the other two 
methods. Moreover, Table [2] suggests that running times for the interior point method do not 
depend too heavily on the type of problem being solved. 
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Figure 3: Mean and median running times (seconds) for three methods on ten problem classes. 

6 Conclusions 

By incorporating a particularly efficient linear-system solution in the Newton update, a generic 
interior point method can perform as well as, or better than, special-purpose methods on convex 
resource allocation problems. This addresses two of the questions posed by Patriksson in his sur- 
vey [25]. First, it shows that it is possible to exploit the particular sparsity of resource allocation 
problems within the setting of a general-purpose optimizer. In particular, we note that this special 
structure could be readily identified during the pre-processing phase of a general solver. Second, 
Patriksson notes that the pegging and breakpoint search methods impose assumptions regarding 
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Table 1: Win percentage of IPM over best breakpoint method 
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Table 2: Mean IPM solution times (in milliseconds) by dimension 



problem class 


n = 10 1 


10 2 


10 3 


10 4 


10 5 


lO 6 


lot-sizing 


8.75 


9.46 


20.17 


125.8 


1022 


8045 


powers over simplex 
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11.02 


24.20 


141.0 


1168 


12281 


p-norm over p-ball 


9.07 


10.81 


21.86 


131.1 


1011 


9009 


inventory management 


9.32 


11.28 


28.08 


147.7 


1404 


11756 


reliability investment 


8.23 


9.60 


21.65 


131.7 


991 


10448 


resource renewal 


8.18 


9.88 


24.00 


131.7 


1118 


12082 


p-norm over r-ball 


8.87 


10.76 


20.80 


135.1 


1073 


9235 


powers s.t. powers 


11.55 


13.66 


36.23 


215.1 


1693 


17989 


quartic s.t. linear 


8.69 


13.58 


24.49 


171.2 


1651 


19390 


log-exponential 


8.88 


12.35 


31.83 


233.7 


2807 


38778 



domain, monotonicity or strict convexity of the objective and/or constraint functions (cf., As- 
sumptions A1-A3 in $2j). Interior point methods require little beyond convexity and second-order 
continuous differentiability. By operating in the interior of the intervals they are largely 

unaffected by asymptotes outside those intervals and could potentially even contend with asymp- 
totes at the endpoints. Also, the implementation used here seems to handle a nonlinear inequality 
constraint just as easily as a linear equality. 

Two main directions warrant further investigation. We have provided computational evidence 
of the interior point method's usefulness, but we do not address theoretical issues of convergence or 
complexity. Also, because interior point methods are not easily warm-started, no conclusions are 
drawn here regarding the use of such methods for applications in which the problem data change 
slightly from call to call, as happens in a branch-and-bound algorithm. 
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Figure 4: Extreme and median running times (seconds) for the interior point method. 
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